###########################################################################
###########################################################################
####Non-state Atrocities in Capital Cities - Mechanism 2 (Regressions) ####
###########################################################################
###########################################################################
library(MASS) 
library(pscl)
library(foreign)
library(stargazer)
library(Zelig)
library(ZeligChoice)

# Set working library
setwd("~/Data/Mechanism 2/")
# Read in data for plotting
m.dat.b <- read.dta("mech2analfull.dta")
# Subset out violent campaigns
m.dat.c <- subset(m.dat.b, nvlcamp<2)

#####################
###Baseline models###
#####################
###Run model on all campaigns
##Remove NAs from data
z.dat.1.a <- na.omit(m.dat.b[,c("ccode", "year", "camp", "lagincidentnonstatefull", "year96", "year97", "year98", "year99", "year00", "year01", "year02", "year03")])
#Run model
z.log.1.a <- zelig(camp ~ lagincidentnonstatefull + 
                     year96 + year97+
                     year98 + year99 + year00 + year01 + year02 + 
                     year03,
                   data=z.dat.1.a, model = "logit")
summary(z.log.1.a)
##Calculate predicted probabilities
#Range of values on IV
inc_range <- seq(min(m.dat.b$lagincidentnonstatefull), max(m.dat.b$lagincidentnonstatefull), 1)
#Simulate first differences
out.z.1.a <- setx(z.log.1.a, lagincidentnonstatefull = inc_range)
z.1.a.p <- sim(z.log.1.a, x=out.z.1.a)
#Plot
pdf("baseallcdc.pdf")
ci.plot(z.1.a.p,  xlab = "Number of Insurgent Atrocities", ylab = "Predicted Probability of Civil Disobedience", ci = c(80, 90, 95))
dev.off()

###Run model without primarily violent campaigns
##Remove NAs from data
z.dat.1.b <- na.omit(m.dat.c[,c("ccode", "year", "camp", "lagincidentnonstatefull", "year96", "year97", "year98", "year99", "year00", "year01", "year02", "year03")])
#Run model
z.log.1.b <- zelig(camp ~ lagincidentnonstatefull + 
                     year96 + year97+
                     year98 + year99 + year00 + year01 + year02 + 
                     year03,
                   data=z.dat.1.b, model = "logit")
summary(z.log.1.b)
##Calculate predicted probabilities
#Range of values on IV
inc_range <- seq(min(m.dat.c$lagincidentnonstatefull), max(m.dat.c$lagincidentnonstatefull), 1)
#Simulate first differences
out.z.1.b <- setx(z.log.1.b, lagincidentnonstatefull = inc_range)
z.1.b.p <- sim(z.log.1.b, x=out.z.1.b)
#Plot
pdf("basenvcdc.pdf")
ci.plot(z.1.b.p,  xlab = "Number of Insurgent Atrocities", ylab = "Predicted Probability of Civil Disobedience", ci = c(80, 90, 95))
dev.off()


#################
###Full models###
#################
###Run model on all campaigns
##Remove NAs from data
z.dat.2.a <- na.omit(m.dat.b[,c("ccode", "year", "camp", "lagpolity2", "loglagrgdppc", "lagcivconf", "lagincidentnonstatefull", "year96", "year97", "year98", "year99", "year00", "year01", "year02", "year03")])
#Run model
z.log.2.a <- zelig(camp ~ lagincidentnonstatefull + lagpolity2 + 
                     loglagrgdppc + lagcivconf +
                     year96 + year97+
                     year98 + year99 + year00 + year01 + year02 + 
                     year03,
                   data=z.dat.2.a, model = "logit")
summary(z.log.2.a)
##Calculate predicted probabilities
#Range of values on IV
inc_range <- seq(min(m.dat.b$lagincidentnonstatefull), max(m.dat.b$lagincidentnonstatefull), 1)
#Simulate first differences
out.z.2.a <- setx(z.log.2.a, lagincidentnonstatefull = inc_range)
z.2.a.p <- sim(z.log.2.a, x=out.z.2.a)
#Plot
pdf("fullallcdc.pdf")
ci.plot(z.2.a.p,  xlab = "Number of Insurgent Atrocities", ylab = "Predicted Probability of Civil Disobedience", ci = c(80, 90, 95))
dev.off()

###Run model without primarily violent campaigns
##Remove NAs from data
z.dat.2.b <- na.omit(m.dat.c[,c("ccode", "year", "camp", "lagpolity2", "loglagrgdppc", "lagcivconf", "lagincidentnonstatefull", "year96", "year97", "year98", "year99", "year00", "year01", "year02", "year03")])
#Run model
z.log.2.b <- zelig(camp ~ lagincidentnonstatefull + lagpolity2 + 
                     loglagrgdppc +  lagcivconf +
                     year96 + year97+
                     year98 + year99 + year00 + year01 + year02 + 
                     year03,
                   data=z.dat.2.b, model = "logit")
summary(z.log.2.b)
#Range of values on IV
inc_range <- seq(min(m.dat.c$lagincidentnonstatefull), max(m.dat.c$lagincidentnonstatefull), 1)
#Simulate first differences
out.z.2.b <- setx(z.log.2.b, lagincidentnonstatefull = inc_range)
z.2.b.p <- sim(z.log.2.b, x=out.z.2.b)
#Plot
pdf("fullnvcdc.pdf")
ci.plot(z.2.b.p,  xlab = "Number of Insurgent Atrocities", ylab = "Predicted Probability of Civil Disobedience", ci = c(80, 90, 95))
dev.off()

